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Abstract 



The Gross-Pitaevskii equation assumes that all (identical) bosons of a condensate 
reside in a single one-particle function. Here, we raise the question whether it always 
provides the best mean-field ansatz for condensates, leading to the lowest mean- 
field ground state energy. To this end, we derive a mean-field approach allowing for 
bosons to reside in several different one-particle functions. The number of bosons in 
\ each of these functions is a variational parameter minimizing the energy. The energy 

• and one-particle functions at these optimal numbers can be determined directly. A 

numerical example is presented demonstrating that the mean-field energy of trapped 
bosons can be below that provided by the Gross-Pitaevskii equation. Implications 



£f) ' are discussed. 
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Non-linear Schrodinger equations are widely used in many areas of physics 
and chemistry. In fact, the mean-field equation of any many-body system is 
a non-linear Schrodinger equation. Usually, however, one speaks of the non- 
linear Schrodinger equation if this equation takes on a particular form. This 
equation is used to describe many phenomena like solitons [1], optical pulses 
in fibres [2], Bose-condensed photons [3] and others. Motivated by famous 
experiments [4,5,6], extensive use of this equation has been made in resent 
years to describe dilute gas Bose-Einstein condensates [7]. In the context of 
such condensates this equation is also referred to as Gross-Pitaevskii (GP) 
equation [8,9]. 

The non-linear Schrodinger equation can be derived as the mean-field equa- 
tion of a system of N identical bosons which interact via a 5-function potential 
W(fi — r 2 ) = \o5(fi — r 2 ), where r, is the position of the i-th particle. It is 
equivalent to the GP equation if we identify the appearing non-linear parame- 
ter with the s-wave scattering length of the particles [7] . At zero temperature 
all bosons can occupy the same spatial one-particle function (orbital) if their 
mutual interaction is neglected. The mean-field description of the interacting 
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system is obtained by assuming the ground state wave function ^ to be a prod- 
uct of identical spatial orbitals ip: \&(r*i, r 2 , . . . , rjy) = ip(fi)ip(f 2 ) • ■ • The 
energy = < > defined as the expectation value of the widely used 

[7,10] Hamiltonian H = h + W, reads 

E GP = N{J ^h^dr + ^J M 4 df},A = Ao(iV-l) (1) 



where h(f) = T + V{r) is the unperturbed one-particle Hamiltonian consisting 
of the kinetic operator T and the external potential V(f). By minimizing this 
energy one readily obtains the GP equation for the optimal orbital ip: 

{h(f) + X\ v (f)\ 2 } v (f) = /icp^f) (2) 



Solution of this eigenvalue equation determines this ip and hop which is called 
the chemical potential [7]. 

The GP equation and energy are very appealing and have been successful in 
explaining experimental observations. We raise the question, however, whether 
the GP approach always provides the best mean-field description of the system. 
We shall demonstrate below that other mean-fields may provide a lower energy 
than (1) does, at least in some cases. The general wave function of N non- 
interacting bosons is a product piifi) ip 2 (r 2 ) • • • Vjv(^v) of orbitals which can 
all be different. Since the bosons are identical, this product must, of course, 
be symmetrized and several bosons may occupy the same spatial orbital. We 
may put n\ bosons in orbital ipi, n 2 in orbital (p 2 and so on. By definition 
the most general mean- field energy is E =< ty\H\ty > where \P is the above 
symmetrized product. For transparency of presentation we confine ourselves 
in the following to two orbitals ipi and ip 2 with particle occupations ri\ and 
n 2 , respectively, where rii +n 2 = N. The extension to more orbitals is trivial. 
The energy now reads 

nArii — 1) f , l4 
E = nihu + n 2 h 22 + A J |y?i| dr + 

A J \<p 2 \ 4 dr + 2A n!n 2 J |(/9i| 2 |(/9 2 | 2 (if (3) 

where ha = J •p*h , pidr is the usual one-particle energy. Minimizing this en- 
ergy with respect to <fi and <p 2 under the constraints that they are orthogonal, 
< y?i|y?2 >= 0, and normalized to unity, < y^ly^ >= 1, we obtain the following 
equations 



{ h(r) + A (ni - l)|v?i(r)| 2 + 2A n 2 |v9 2 (f)| 2 } Vl {f) 
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= /in <Pi(r) +//12 ^(r) 
{ h(r) + A (n 2 - l)|v? 2 (r)| 2 + 2\ n 1 \<p 1 (f)\ 2 } cp 2 (r) = 

= H22 ¥2 (r) +^21^1 (r) (4) 

for the optimal orbitals. These equations should not be confused with those for 
condensates made of two types of bosons. The /i^ are the Lagrange parameters 
due to the above mentioned constraints. Of course, the energy E and the 
orbitals y?i,y? 2 depend on the particle occupation ri\. Consequently, we may 
treat n\ as a variational parameter and search for its optimal value which 
minimizes the energy. We shall return to this central point below. 

Here, we briefly make contact with the problem of fragmentation often dis- 
cussed in the literature, see refs. [11,12] and references therein. A condensate 
is fragmented if its reduced one-body density matrix has two or more macro- 
scopic eigenvalues. Until now no fragmentation has been found for confined 
condensates, i.e., for condensates in an external potential. Moreover, it has 
been shown that fragmentation cannot take place in a harmonic external po- 
tentials [13]. Our ansatz discussed above in eqs.(3,4) can describe fragmenta- 
tion on the mean-field level. If the particle occupations n\ and n 2 minimizing 
the energy are both macroscopic for large N, fragmentation takes place. 

Let us analyze eqs.(3,4). Using eq.(l) it follows that the energy needed to 
remove a particle without changing the orbital tp is, within the GP mean- 
field, given by Eqp(N) — E GP (N — 1) = /igp- Analogously, we can compute 
the energy needed to remove a particle from the orbital y?i and that from 
the orbital (p 2 in the framework of the present mean-field. Using eqs.(3,4) and 
remembering that cpi and f 2 are orthogonal, we readily find: /in = E(ni,n 2 ) — 
E(ni — 1, n 2 ) and /i 22 = E{n\, n 2 ) — E(n\, n 2 — 1). In other words, /in and /i 22 
can be viewed as chemical potentials of the {y?i}- and {<y9 2 }-particle manifolds, 
respectively. 

Another interesting property of (3,4) is that, unless n\ = n 2 = N/2, these 
equations are not invariant to rotations (unitary transformations) of the pair 
(ipi, <p 2 ). This is due to the fact that ipi and <p 2 appear a different number of 
times in the wave function in contrast to the situation for identical fermions 
where the mean-field equations (Hartree-Fock) are indeed invariant to rota- 
tions [14]. As an important consequence the off-diagonal Lagrange parameters 
/ii 2 and fi 2 i can not, in general, be removed from eq.(4) by an unitary trans- 
formation. Finally we mention that ni/i i2 = n 2 /i 2 i. 

The factor 2A nin 2 appearing in the energy (3) and the corresponding term 
in (4) deserve attention. If the bosons residing in </?i would be distinguishable 
from those living in tp 2 , this factor will reduce to Aonin 2 . The remaining factor 
Ao^i^ 2 appears because of the exchange interaction between the two manifolds 
of identical bosons. Concentrating on many bosons, i.e. ni,n 2 ^> 1, equations 
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(3,4) become more transparent. Then we may introduce the orbitals <pi = 
(ni/N) 1 / 2 (fi and (fi 2 — (^2 /N) 1 ^ 2 ip 2 which are, of course orthogonal, but 
now normalized to rii/N and n 2 /N, i.e., to the fraction of particles residing 
in tpi and ip 2 , respectively. The energy (3) can now be expressed as 

E = E c i + E ex 

A f 

E d = h u + h 22 + -J p 2 (f)dr (5) 
E ex = \f |0i| 2 |0 2 | 2 cF 



where hu is defined as hu in (3), but with the orbitals (pi and 

p = Mv?i| 2 + n 2 \<p 2 \ 2 } /N = |0!| 2 + |0 2 | 2 

is the density per particle. It is easy to see that the 'classical' energy E d is 
invariant to unitary transformation of the pair (<pi,(p 2 ) while the exchange 
energy E ex is not. It is the latter energy which makes the situation interest- 
ing. Without the exchange energy also the equations for the optimal orbitals, 
eq.(4), would be invariant to unitary transformations. 

In the following we demonstrate that the mean-field energy (3) can be superior 
to the GP energy (1). Our example consists of an one-dimensional attractive 
condensate (A < 0) in an external symmetric double-well potential. We choose 
as potential (a > \x \) 
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V(x) = c + -lux — ua x — x + a (6) 



where the bottoms of the well are at x — ±xq- The parameter a is chosen such 
that V(±x ) = 0. The corresponding kinetic energy reads T = — f ^2 implying 
that a; is a dimensionless coordinate. Since the potential is symmetric, the 
ground state of h is a gerade function and the first excited state is an ungerade 
function of x. We, therefore, consider in our mean-field tpi and ip 2 to be gerade 
and ungerade functions of x, respectively. As can be seen from eq.(4) one finds 
M12 = M21 = because tpi and tp 2 possess different spatial symmetries. The 
equations for the optimal orbitals have been evaluated numerically, by solving 
the resulting eigenvalue equations self-consistently using DVR [15] and finite- 
difference [16] methods. Both methods have led to identical results. 

Using the same numerical procedures we have also solved the GP equation (2) 
and determined ip. Having found the optimal orbitals, we have computed the 
mean-field energies E and Eqp- Figure 1 shows the difference of the energies 
per particle (E — Egp)/N as a function of the fractional occupation n\/N 
for several values of the coupling constant A = A iV (note that N ^> 1). At 
ni/N — 1 the energy E, of course, coincides with the GP energy Egp for all 
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values of A. For A = 0.5 the energy E grows as n\/N is decreased. This is 
also true for A = —0.5, but the growth of E is less pronounced. Until about 
A ps —1.0, E > Eqp implying that the GP mean field is the best mean field 
and no particle occupy (p 2 . At A < —1.0 the scenario changes as the energy 
E drops below E G p as a function of the fractional occupation n\/N for a 
given value of A. The optimal mean field is obtained at the minimum of the 
curve E = E{n\/N). This minimum is about n\/N = 3/4 for A = —1.5 and 
approaches n-y/N = 1/2 continuously as A is further decreased to A = —4.0. 
The energy gain per particle due to the present mean field is substantial; for 
the same potential the energy difference between the ground and first excited 
states of h = T + V is ~ 0.18 and the height of the barrier separating the two 
wells is ~ 1.0. 

The results are better understood by inspecting fig. 2 which shows the density 
per particle p(x) for A = —2.5 and rii/N = 0.6. As can be seen in fig.l, this 
value of ni/N corresponds to the minimum of the energy for A = —2.5. The 
GP density per particle |<^(x)| 2 is also shown in fig. 2 for comparison. It is 
eye catching that much of |v?(a;)| 2 has moved from the center of the potential 
x 0.0 to its wells as an impact of the present mean field. Furthermore, the 
two peaks of the density have narrowed. The density p{x) has substantially 
localized in comparison to the GP density. Since the condensate is attractive 
(A < 0), localization is well-known to lead to a lowering of the energy [10]. 
What is surprising is that localization beyond the symmetry preserved GP 
density can be achieved in the framework of a mean-field approach without 
breaking the symmetry and without the necessity to go beyond the mean- 
field world. We remark that the system should have more than two bosons to 
achieve this localization: because of the different spatial symmetries of ipi and 
ip 2 , their product is not gerade as the ground state is expected to be. If more 
bosons are present, an even number of them can reside in ip 2 . 

The energy E has been discussed above as a function of rt\. The particle oc- 
cupation ri\ is a variational parameter and we may ask whether the minimum 
of E = E{n\) can be determined directly without computing first the energy 
as a function of n\ and subsequently searching for the minimum. The answer 
is positive. We define a new function 

r(n 1 ,n 2 ) = E(n 1 ,n 2 ) + v(ni + n 2 - N) (7) 

of rii and n 2 and the constraint n\ + n 2 = N is represented by the Lagrange 
parameter v. At its extrema dr/drii = dr/dn 2 = 0, leading to two relations, 
one of which fixes v and other gives rise to the equality dE/dni = dE/dn 2 . 
Taking the derivative of E in (3) with respect to n\ and n 2 for n\,n 2 ^> 
1 recovers the expressions for /in and p 22 determined from (4), i.e., pa = 
dE/drii, i — 1, 2. Consequently we find that the two quantities pn and p 22 are 
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equal at the extremum of E: 



/in = H22 = H (8) 



This is a relevant finding for understanding the concept of a condensate and, 
furthermore, allows us to transform the basic eqs.(4) at the optimal occupation 
n\ as discussed in the following. 

We now introduce the new orbitals 



^i(r) = (^) 1/ Vi(r) + (|) 1/ V 2 (r) 

^ 2 (r) = (|) 1/ V 2 (r)-(^) 1/ Vi(r) (9) 

and express (4) by them. Confining ourselves to the minimum of E and using 
the basic relation (8), the resulting equations for tpi and tp 2 simplify consider- 
ably: 



{ h(f) + ^A|^| 2 + ^MH 2 } V»i= (a* + P) & 

{ h(f) + ^A|^ 2 | 2 + \x\H 2 } ^2 = (jm - p) ife. (10) 



The equations (4) for the orbitals </?i and ip 2 are, in general, not eigenvalue 
equations. The corresponding equations for ipi and ijj 2 at the optimum occupa- 
tion n\ are eigenvalue equations. The eigenvalues are (/i+/t) and (/i— /t), where 
/t is the condensate's chemical potential (see (8)) and /t = (ni/n2) 1//2 /ii2 = 
(n 2 /ni) 1/2 /i2i- 

At first sight it is intriguing to see that the new eqs.(10) do not exhibit the 
optimal rii and n 2 . Nevertheless, these occupations can be deduced from these 
equations. To this end it is relevant to notice that (9) is not a unitary transfor- 
mation. Indeed, ipi and ip2 are normalized to unity, < >=< ip 2 1^2 >= 1) 
but they are not orthogonal. Their scalar product is 

< >= (id 



and determines the optimal occupations. One can solve the two eigenvalue 
equations (10) self-consistently and compute the optimal n\ and n 2 from the 
final result. We have tested this procedure for several numerical examples in- 
cluding the above discussed double-well case. One directly obtains the minima 
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of the energy curves shown in fig. 1. A practical hint: numerical stability is en- 
hanced if at each iteration step use is made of the relationship (9) between ip 1 
and ip 2 - 

Let us conclude. By expressing the wave functions as a symmetrized product of 
single-particle functions (orbitals), equations for the mean-field energy and the 
optimal orbitals are derived. Several bosons can reside in each of these orbitals 
ifi and it is a basic property of the present approach that these occupation 
numbers are variational parameters which can be used to minimize the 
energy. An explicit numerical example is presented demonstrating that the 
mean-field energy using two orbitals can be lower than the GP energy. Since 
in the absence of interaction between the particles, all particles reside in a 
single orbital to form the ground state of the condensate, it is clear that the 
present approach reproduces the GP result for weak interactions. Beyond some 
critical interaction strength A = A N obtained, for instance, by increasing the 
number N of particles in the condensate, the mean-field energy can drop below 
the GP one. This critical value depends strongly on the external potential. By 
varying the shape and parameters of this potential one can lower or raise the 
interesting range of A substantially. 

In our simple example it has been shown that the mean-field energy can drop 
below the GP energy for a continuous range of the occupation rii. The system 
then consists of two (or more) subsystems each possessing its own chemical 
potential. This may contradict the picture one usually has of a condensate. It 
is, therefore, very interesting to note that these, in general different, chemical 
potentials become identical at the optimal occupations rii which minimize the 
energy E({rii}) restoring thereby the picture of a condensate. This finding 
leads to another highlight of this work: self-consistent eigenvalue equations for 
the optimal orbitals at the optimal occupations could be derived. 

It should be kept in mind that the occupation numbers are variational pa- 
rameters in the present mean-field and, therefore, the solution of eq.(10) will 
tell us via eq.(ll) which the optimal occupation numbers are. In those cases 
where one finds that rii/N = 1, then all bosons reside in a single orbital, our 
ansatz reproduces the GP result and the GP mean field is thus the best mean 
field. The GP mean field is an excellent ansatz reflecting our physical expecta- 
tion that in the ground state all bosons would like to occupy the same orbital 
in order to minimize the energy. It is, therefore, surprising that there exist 
counter examples at all. Here, we have shown that an attractive condensate in 
an one-dimensional symmetric double-well potential provides such a counter 
example. This finding is not restricted to the specific shape of the external 
symmetric double-well potential used here. Other counter examples can be 
expected. For instance, an attractive condensate in an one- dimensional sym- 
metric multi-well potential. Here, however, as many orbitals must be included 
in the mean-field ansatz as there are wells. That macroscopic condensation can 
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be described by one-dimensional equations has been discussed in the litera- 
ture, see, e.g. [10] and references therein. Clearly, it would be very interesting 
to study condensates in higher dimensions. Of course, the same equations as 
derived here hold also in this case, but their full numerical evaluation is more 
involved and such a study is left to the future. We note that in the present 
numerical example it is the symmetry of the problem which fails the GP mean 
field. In other examples the reason that GP is not the best mean field may 
be of different origin. We suspect, however, that if one finds rii/N to deviate 
from 1, it is a helpful warning that one might have to go beyond the mean-field 
approach as such in order to describe the underlying physics properly. 

Being non-linear equations, eqs. (4) and (10) may have several solutions. We 
have encountered examples where the lowest-energy solution is characterized 
by rii/N = 1, but other low- lying solutions exist with n-y/N deviating from 
I. The first solution, of course, corresponds to the system's ground state. In 
the latter solutions the bosons reside in two orbitals and their existence could 
be an indication for collective excited self- consistent states of the condensate. 
We remind that according to eq. (8) there exists only a single value of the 
chemical potential for each of these excited states making them attractive to 
further investigations. 

Finally, we would like to mention very briefly that (a) time-dependent mean- 
field equations corresponding to the stationary equations discussed here are 
straightforwardly derived, and (b) the present mean-field equations provide 
an excellent starting point of investigating condensates beyond the mean-field 
approach. 

Helpful discussions with Dieter Meyer and Ofir Alon are gratefully acknowl- 
edged. 
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Fig. 1. The difference between the present mean-field energy per particle E(n\)/N 
and the GP energy per particle Eqp/N is shown for several values of A = Ao N 
as a function of the fractional occupation n±/N. It is seen that the energy E(n±) 
may drop substantially below the GP energy. At n\/N = 1 both energies, of course, 
coincide. The external potential (4) has been used with xq = 1.5 and a = 1.5015. A 
and the energies are given in units of u> 
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Fig. 2. The density per particle p = (ni|</?i| 2 + ri2\ip2\ 2 )/N is shown for one of the 
examples of fig.l (A = —2.5 and n\/N = 0.6) and compared with the corresponding 
GP density per particle pep = M 2 - 
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